"\n"
"// This function for analytical evaluation of particle density integral is \n"
"// provided by Eric Bruneton\n"
"// http://www-evasion.inrialpes.fr/Membres/Eric.Bruneton/\n"
"//\n"
"// optical depth for ray (r,mu) of length d, using analytic formula\n"
"// (mu=cos(view zenith angle)), intersections with ground ignored\n"
"float2 GetDensityIntegralAnalytic(float r, float mu, float d, float EarthRadius, float4 ParticleScaleHeight) \n"
"{\n"
"    float2 f2A     = sqrt(ParticleScaleHeight.zw * 0.5 * r);\n"
"    float4 f4A01   = f2A.xxyy * float2(mu, mu + d / r).xyxy;\n"
"    float4 f4A01s  = sign(f4A01);\n"
"    float4 f4A01sq = f4A01 * f4A01;\n"
"    \n"
"    float2 f2X;\n"
"    f2X.x = f4A01s.y > f4A01s.x ? exp(f4A01sq.x) : 0.0;\n"
"    f2X.y = f4A01s.w > f4A01s.z ? exp(f4A01sq.z) : 0.0;\n"
"    \n"
"    float4 f4Y =\n"
"        f4A01s \n"
"        / (2.3193 * abs(f4A01) + sqrt(1.52 * f4A01sq + 4.0)) \n"
"        * float3(1.0, exp(-ParticleScaleHeight.zw * d * (d / (2.0 * r) + mu)) ).xyxz;\n"
"\n"
"    return \n"
"        sqrt((6.2831 * ParticleScaleHeight.xy) * r) \n"
"        * exp((EarthRadius - r) * ParticleScaleHeight.zw) \n"
"        * (f2X + float2(f4Y.x - f4Y.y, f4Y.z - f4Y.w));\n"
"}\n"
"\n"
"\n"
"float3 GetExtinctionUnverified(float3 f3StartPos,\n"
"                               float3 f3EndPos,\n"
"                               float3 f3EyeDir,\n"
"                               float3 f3EarthCentre,\n"
"                               float  fEarthRadius,\n"
"                               float4 f4ParticleScaleHeight)\n"
"{\n"
"#if 0\n"
"    float2 f2ParticleDensity = IntegrateParticleDensity(f3StartPos, f3EndPos, f3EarthCentre, 20);\n"
"#else\n"
"    float r = length(f3StartPos-f3EarthCentre);\n"
"    float fCosZenithAngle = dot(f3StartPos-f3EarthCentre, f3EyeDir) / r;\n"
"    float2 f2ParticleDensity = GetDensityIntegralAnalytic(r, fCosZenithAngle, length(f3StartPos - f3EndPos), fEarthRadius, f4ParticleScaleHeight);\n"
"#endif\n"
"\n"
"    // Get optical depth\n"
"    float3 f3TotalRlghOpticalDepth = g_MediaParams.f4RayleighExtinctionCoeff.rgb * f2ParticleDensity.x;\n"
"    float3 f3TotalMieOpticalDepth  = g_MediaParams.f4MieExtinctionCoeff.rgb      * f2ParticleDensity.y;\n"
"        \n"
"    // Compute extinction\n"
"    float3 f3Extinction = exp( -(f3TotalRlghOpticalDepth + f3TotalMieOpticalDepth) );\n"
"    return f3Extinction;\n"
"}\n"
"\n"
"float3 GetExtinction(float3 f3StartPos,\n"
"                     float3 f3EndPos,\n"
"                     float3 f3EarthCentre,\n"
"                     float  fEarthRadius,\n"
"                     float  fAtmTopRadius,\n"
"                     float4 f4ParticleScaleHeight)\n"
"{\n"
"    float3 f3EyeDir = f3EndPos - f3StartPos;\n"
"    float fRayLength = length(f3EyeDir);\n"
"    f3EyeDir /= fRayLength;\n"
"\n"
"    // Compute intersections of the view ray with the atmosphere and the Earth\n"
"    float4 f4Isecs;\n"
"    GetRaySphereIntersection2( f3StartPos, f3EyeDir, f3EarthCentre,\n"
"                               float2(fEarthRadius, fAtmTopRadius), \n"
"                               f4Isecs);\n"
"    float2 f2RayEarthIsecs  = f4Isecs.xy;\n"
"    float2 f2RayAtmTopIsecs = f4Isecs.zw;\n"
"\n"
"    // If the ray misses the atmosphere, there is no extinction\n"
"    if (f2RayAtmTopIsecs.y < 0.0)\n"
"    {\n"
"        return float3(1.0, 1.0, 1.0);\n"
"    }\n"
"\n"
"    // Limit the ray length by the distance to the top of the atmosphere\n"
"    fRayLength = min(fRayLength, f2RayAtmTopIsecs.y);\n"
"\n"
"    if (f2RayEarthIsecs.x > 0.0)\n"
"    {\n"
"        // The ray hits the Earth surface\n"
"        fRayLength = min(fRayLength, f2RayEarthIsecs.x);\n"
"    }\n"
"\n"
"    // Update the end point\n"
"    f3EndPos = f3StartPos + f3EyeDir * fRayLength;\n"
"\n"
"    // If the ray starts outside of the atmosphere, move the starting point to the intersection\n"
"    f3StartPos += f3EyeDir * max(f2RayAtmTopIsecs.x, 0.0);\n"
"\n"
"    return GetExtinctionUnverified(f3StartPos, f3EndPos, f3EyeDir, f3EarthCentre, fEarthRadius, f4ParticleScaleHeight);\n"
"}\n"
